#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <time.h>
#include <pthread.h>

#define ITER 30  //Numero de iteracoes que serao executadas

//Declaracao das variaveis utilizadas durante as iteracoes. Declaradas publicamente para que as threads tenham acesso.
mpf_t a[ITER]; 
mpf_t b[ITER];
mpf_t t[ITER];
mpf_t p[ITER];
pthread_mutex_t fA[ITER];
pthread_mutex_t fB[ITER];
mpf_t pi;

FILE *piFile;

//Funcoes que calculam os valores de cada variavel a cada iteracao
void *f_calculaA(void *empty) 
{ 
	long i;
	for (i = 1 ; i < ITER ; i++)
	{
		pthread_mutex_lock(&fB[i-1]);
		pthread_mutex_unlock(&fB[i-1]);
		mpf_add(a[i],a[i-1],b[i-1]);
		mpf_div_ui(a[i],a[i],2); //an+1 = (an + bn)/2
		pthread_mutex_unlock(&fA[i]);
	}
}

void *f_calculaB(void *empty)
{
	long i;
	for (i = 1 ; i < ITER ; i++)
	{
		pthread_mutex_lock(&fA[i-1]);
		pthread_mutex_unlock(&fA[i-1]);
		mpf_mul(b[i],a[i-1],b[i-1]);
		mpf_sqrt(b[i],b[i]);
		pthread_mutex_unlock(&fB[i]);
	}
}

void *f_calculaT(void *empty)
{
	long i;
	for (i = 1 ; i < ITER ; i++)
	{
		pthread_mutex_lock(&fA[i]);
		pthread_mutex_unlock(&fA[i]);
		mpf_sub(t[i],a[i-1],a[i]);
		mpf_mul(t[i],t[i],t[i]);
		mpf_mul(t[i],t[i],p[i-1]);     
		mpf_sub(t[i],t[i-1],t[i]);  // tn+1 = tn - pn*(an-an+1)^2
	}
}

int main (int argc, char *argv[]) 
{
	//Declaracao de cada uma das threads que serao executadas
	pthread_t calculaA;
	pthread_t calculaB;
	pthread_t calculaT;

	//Abertura do arquivo onde sera armazenado o valor final calculado
	piFile = fopen("out_gaussLegendrePar.txt","w+");

	if(piFile == NULL) 
	{
		return EXIT_FAILURE;
	}

	//Variaveis usadas para calcular tempo levado pelo algoritmo
	time_t begin, end;
	double time_spent;
	long i;

	printf("Initializing variables...");
	printf("\n");
	
	//Setar precisão (usar todas as 10 milhões de casas) e inicializar variaveis do metodo
	mpf_set_default_prec(40000000);
	
	mpf_init_set_ui(a[0],1);
	
	mpf_init(b[0]);
	mpf_sqrt_ui(b[0],2);
	mpf_ui_div(b[0],1,b[0]);
	
	mpf_init_set_str(t[0],"0.25",10);
	
	mpf_init_set_ui(p[0],1);
	
	mpf_init(pi);
	
	for (i = 1 ; i < ITER ; i++)
	{
		mpf_init_set_si(a[i],-1);
		mpf_init_set_si(b[i],-1);
		mpf_init_set_si(t[i],-1);
		mpf_init_set_si(p[i],-1);
		mpf_mul_ui(p[i],p[i-1],2);
		pthread_mutex_init(&fA[i],NULL);
		pthread_mutex_lock(&fA[i]);
		pthread_mutex_init(&fB[i],NULL);
		pthread_mutex_lock(&fB[i]);
	}
	
		
	
	printf("Starting...");
	printf("\n");
	time(&begin);

	pthread_create(&calculaA,NULL,f_calculaA,NULL);
	pthread_create(&calculaB,NULL,f_calculaB,NULL);
	pthread_create(&calculaT,NULL,f_calculaT,NULL);
	
	pthread_join(calculaA,NULL);
	pthread_join(calculaB,NULL);
	pthread_join(calculaT,NULL);
	
	//Calcula valor de pi depois das iteracoes
	mpf_add(pi,a[ITER-1],b[ITER-1]);
	mpf_mul(pi,pi,pi);
	mpf_div_ui(pi,pi,4);
	mpf_div(pi,pi,t[ITER-1]); //pi = ((an + bn)^2)/(4*tn)

	//Finaliza o cronometro e imprime o tempo de execucao
	time(&end);
	time_spent = difftime(end, begin);
	printf("Time: %lf seconds - Iterations: %ld\n",time_spent, i);
	
	//Escreve no arquivo o valor calculado de pi
	mpf_out_str(piFile,10,250000,pi);

	//Fecha o arquivo e limpa as variaveis
	fclose(piFile);
	//mpf_clears(pi,aN,bN,tN,pN,aNP1,bNP1,tNP1,pNP1,NULL);
}
